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Abstract 

The self-organized Monte Carlo simulations of 2D Ising 
ferromagnet on the square lattice are performed. The 
essence of devised simulation method is the artificial 
dynamics consisting of the single-spin-flip algorithm of 
Metropolis supplemented by the random walk in the 
temperature space. The walk is biased to the criti- 
cal region through the feedback equation utilizing the 
memory-based filtering recursion instantly estimating 
the energy cumulants. The simulations establish that 
the peak of the temperature probability density func- 
tion is located nearly the pseudocritical temperature 
pertaining to canonical equilibrium. In order to elimi- 
nate the finite-size effects, the self-organized approach 
is extended to multi-lattice systems, where feedback is 
constructed from the pairs of the instantaneous run- 
ning fourth-order cumulants of the magnetization. The 
replica-based simulations indicate that several prop- 
erly chosen steady statistical distributions of the self- 
organized Monte Carlo systems resemble characteristics 
of the standard self-organized critical systems. 
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1 Introduction 

The Monte Carlo (MC) simulation methods are nonper- 
turbative tools of the statistical physics developed hand 
in hand with increasing power of nowadays computers. 
The benchmark for testing of MC algorithms represents 
the exactly solvable Ising spin model [1]. Between al- 
gorithms applied to different variants of this model, the 
using of the single-spin- flip Metropolis algorithm [2, 3] 
prevails due to its simplicity. Nevertheless the principal 
problems have emerged as a consequence of the accu- 
racy and efficiency demands especially for the critical 
region. 

As a consequence of this several methods enhancing 
MC efficiency have been considered. The procedures of 
great significance are: finite-size-scaling relations and 
renormalization group based algorithms [4] - [7], clus- 
ter algorithms [8, 9] lowering the critical slowing down, 
histogram and reweighting techniques [10] interpolat- 
ing the stochastic data, and multicanonical ensemble 
methods [11] overcoming the tunneling between coex- 
isting phases at 1st order transitions. Even with the 
mentioned modifications and related MC versions, a la- 
borious and human assisted work is needed until a sat- 
isfactory accuracy of results is achieved. This is trivial 



reason why utilization of the self-organization principles 
have attracted recent attention of MC community. 

In the present paper we deal with the combining 
of the self-organization principles and MC dynamics. 
Of course, this effort has a computational rather than 
physical impact. Our former aim was to design the 
temperature scan seeking for the position of the critical 
points of the lattice spin systems. But this original aim 
has been later affected by the general empirical idea of 
the self-organized criticality (SOC) [12], originally pro- 
posed as a unifying theoretical framework describing 
a vast class of systems evolving spontaneously to the 
critical state. The SOC examples are sand pile, forest- 
fire [13] and game-of-life [14] models. The SOC prop- 
erty can be defined through the scale-invariance of the 
steady state asymptotics reached by the probability den- 
sity functions (pdf 's) constructed for spatial and tem- 
poral measures of dissipation events called avalanches. 
The dynamics of standard SOC systems is governed by 
the specific noncquilibrium critical exponents linked by 
the scaling relations [15] in analogy to standard phase 
transitions. 

Notice that dynamical rules of standard SOC systems 
are functions of the microscopic parameters uncoupled 
to the global control parameters. On the contrary, 
in the spin systems governed by the single-spin-flip 
Metropolis dynamics, the spins are flipped according to 
prescriptions depending on their neighbors, but also on 
the global or external control parameters, like temper- 
ature or selected parameters of the Hamiltonian. The 
modification, by which the MC spin dynamics should 
be affected to mimic the SOC properties, is discussed 
in [16]. In agreement with [17], any such modification 
needs the support of nonlinear feedback mechanism en- 
suring the critical steady stochastic state. It is clear 
that the feedback should be defined in terms of MC in- 
stant estimates of the statistical averages admitting the 
definition of critical state. 

The real feedback model called probability- changing 
cluster algorithm [22] was appeared without any ref- 
erence to the general SOC paradigm. The alternative 
model was presented in [18], where temperature of the 
ferromagnet Ising spin system was driven according to 
recursive formula corresponding to general statistical 
theory [19]. This example was based on the mean mag- 
netization leading to series of the temperature moves 
approaching the magnetic transition. Despite the suc- 
cess in the optimized localization of critical temperature 
of the Ising ferromagnet, the using of term SOC seems 
to be not adequate for this case. The reason is the ab- 



sence of the analysis of spatio-temporal aspects of MC 
dynamics, which can be considered as a noncanonical 
equilibrium [20, 21], due to residual autocorrelation of 
the sequential MC sweeps. 

Regarding the above mentioned approaches, the prin- 
cipal question arises, if the MC supplemented by feed- 
back, resembles or really pertains to branch of the stan- 
dard SOC models which are well-known from the Bak's 
original paper [12] and related works. 

The plan of our paper is the following. Sec. 2 is in- 
tended to generalization of the averaging relevant for 
the implementation of the self-organization principles. 
The details of feedback construction based on the tem- 
perature gradient of specific heat are discussed in Sec. 3. 
These proposals are supplemented by the simulation re- 
sults carried out for 2D Ising ferromagnet. The details 
of the multi-lattice self-organized MC simulations, sta- 
bilizing true critical temperature via fourth-order mag- 
netization cumulants, are presented in Sec. 4. The over- 
simplified mean-field model of self-organized MC algo- 
rithm is discussed in Sec. 5. Several universal aspects of 
the self-organized MC dynamics are outlined by replica 
simulations in Sec. 6. Finally, the conclusions are pre- 
sented. 



2 The running averages 

As we already mentioned in introduction, the mecha- 
nism by which many body system is attracted to the 
critical (or pseudocritical) point should be mediated by 
the feedback depending on the instantaneous estimates 
of the statistical averages. In this section we introduce 
the running averages important to construct the proper 
feedback rules. 

Consider the MC simulation generating the sequence 
of configurations {X t i, t' = 1,2,..., t} according to 
importance sampling update prescription of Metropo- 
lis [3] producing the canonical equilibrium Boltzmann 
distribution as a function of the constant temperature 
T. The estimate of the canonical average (A) t of some 
quantity A 
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can simply calculated from the series of t sampled real 
values Af = A(X t >) reweighted by w t ensuring the triv- 
ial normalization 
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The summation given by Eq.(l) is equivalent to recur- 
rence 
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(1 - w t )(A)t-i,T + w t A t 
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showing how the average changes due to terminal con- 
tribution A t . Consider the generalized averaging, where 
w t is replaced by the constant parameter < r\ <C 1 
which is independent of t: 
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where (1 — 77)' <C 1 if the time of averaging is sufficiently 
large r\ <C 1 <C t. It should be remarked that the gen- 
eralized averages labeled by (A), n .t,T are equivalent to 
gamma filtered [23] fluctuating inputs A t . Note that the 
term gamma originated from the analytic form of the 
weight w Vi t,f- The filtering represents the application 
of the selection principle suppressing the information 
older than the memory depth oc I/77. 

3 The specific heat feedback 

To attain the critical region self-adaptively we construct 
the temperature dependent feedback changing the tem- 
perature in a way enhancing extremal fluctuations. The 
running estimates of averages are necessary to predict 
(with share of the uncertainty) the actual system po- 
sition and the course in the phase diagram leading to 
critical point. 

The pseudocritical temperature T C (N) of some finite 
system consisting of N degrees of freedom is defined by 
the maximum C(N, T C (N)) of the specific heat C(N, T). 
To form an attractor nearly T C (N), we propose the fol- 
lowing dynamics of the temperature random walker 



T t +r t A sign(ff) 



biased by the gradient 
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(A) v ,t.T = (1 - v)(A) v ,t-i,T + f]A t . 
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where C is MC estimate of the specific heat. From 
that follows that the energy fluctuations are controlled 
by the temperature representing the additional slowly 
varying degree of freedom. It is assumed here and in 
further that T t remains constant during N random mi- 
croscopic moves (IMC step per N). The sign function 
occurring in Eq.(9) is used to suppress the extremal 
fluctuations of Ff causing the unstable boundless be- 
havior of T t . From several preliminary simulations it 
can be concluded that the replacement of sign function 
by some smooth differcntiable function (e.g. tangent 
hyperbolic or arcus tangent) seems to be irrelevant for 
keeping of smaller dispersion of T t . The non-constant 
temperature steps are constrained by \T t +\ — T t \ < A 



due to action of the pseudorandom numbers r t drawn 
from the uniform distribution within the interval (0, 1). 

Very important for further purposes is the quasi- 
equilibrium approximation (A) t ~ (-A)n t r t justified un- 
der the restrictions t cq ,A f] -C 1, f cq iA <C 1, where 
t e q,A is the equilibration time of A. With help of this 
approximation, the running averages from Eq.(4) can 
be generalized using 

(A) v ,t.T t = (1 - v)(A)r h t-l,T t - 1 + W t A t>Tt , (11) 

which allow the averaging under the slowly varying tem- 
perature. Here T t is the temperature for which the last 
sample A tt T t — A(X t ) is calculated. 

For later purposes we also introduce the zero passage 
time ta pertaining to A. The time is defined as a mea- 
sure of the stage during which the sign(A t ) is invariant. 
More formal definition of ta requires the introducing of 
two auxiliary times t', t". The first time t' defines the 
instant, where 



A t >-iA t , <0, 
whereas t" counts the events for which 

Af+t»-i Af+t» > , t" = 1,2,... ta > 1 • 
The counting is finished for t" = ta if 



t'+TA Af+TA + 1 



< 0. 
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To be thorough, the conditions A t i_\A t i < 0, 
A t > A t '+i < should be assumed for definition of 
ta = 1- Thus, using Eqs.(12)-(14) the original sequence 
{A t }t=i can be transformed to the sequence of passage 

times {r^ fc) }fe=i. 

From the above definition it follows that random walk 
in temperature is unidirectional (the sign of T t+ \ — T t 
is ensured) for t" = 1,2, ...,Tpc. The arithmetic 
average of t f c is related to temperature dispersion 
((ST) 2 ) = (A 2 /3)(t f c), where (t f c) is calculated from 

{T F k c}k=i- The standard formula providing C(T) is the 
fluctuation-dissipation theorem (in fee units) 
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Here the specific heat is expressed in terms of the en- 
ergy E cumulants (E) t , (E 2 )t- In the frame of the 
quasi-static approximation it is assumed: T ~ T tl 
(E) t .T ~ (£%,t,T t , (E 2 ) t ,T ^ (E 2 ) v , t ,T t - Subsequently, 
using properties of energy cumulants with equilibrium 
Boltzmann weights, the temperature derivative of C can 
be approximated by 



(E 3 ) v .t,T t - 3(E) v ,t,T t (E 2 ) v . 
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The Ising ferromagnet is simulated in further to study 
the effect of feedback defined by Eqs.(9), (10) and (16). 
However, it is worthwhile to note that many of the pre- 
sented results are of general relevance. Given the spin 
system X = {si}f =1 , s, = ±1/2 placed at N — L 2 sites 
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Figure 1: The transient regime of T t obtained for sev- 
eral different initial conditions. Simulated for L = 10, 
i] = 10~ 3 , A = 10~ 4 and identical initial values of cu- 
mulants. 
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Figure 2: The pdf distributions obtained for rj — 10 3 , 
A = 10~ 4 , L = 10: a) the pdf of temperature with the 
peak located nearly T C (N) ~ 0.5868 (dispersion 0.018); 
b) the simulation reveals the flat tails of pdf of Ff . 



of the square L x L lattice with the periodic bound- 
ary conditions, the Ising Hamiltonian can be defined in 
exchange coupling units 



E 



(17) 



where nn means that summation running over the spin 
nearest neighbors. 

In general, the dynamics of SOC systems exhibits 
two distinct regimes. During the transient regime the 
proximity of critical or pseudocritical point is reached. 
The second steady regime is called here the noncanon- 
ical equilibrium in analogy with [21]. In this regime 
the attraction to critical point is affected by the crit- 
ical noise. This general classification is confirmed by 
our results shown in Figs. 1-3. In Fig.l we see the 
stochastic paths of T t pertaining to different initial 
values of energy cumulants and spin configurations. 
The paths are attracted by T C (N) with some uncer- 
tainty in noncanonical equilibrium. For the sufficiently 
narrow steady pdf's, T C (N) can be approximated by 
T C (N) ~ N~ v x Tt, where N av is the number of in- 

puts. The stationary pdf of T t walk is shown in Fig. 2a 
and pdf of Ff with non-Gaussian fiat tails is depicted 
in Fig.2b. 

The alternative quantity capable for the characteri- 
zation of the noncanonical equilibrium is the autocor- 
relation 



— Vf; 
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c v c 
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The simulations results depicted in Fig. 3a evidenced 
that minimum time for which the anticorrelation (K T < 
0) occur is of the order 0((t>c)). As we see from 
Fig. 3b, the power-law dependence can be identified 
within the region of vanishing t f c . More profound dis- 
cussion of this fact is presented in Sec. 6. 



4 Mult i- lattice simulations 
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Figure 3: The simulation for L = 10, r\ = 1CP 3 , A = 
10~ 4 ; a) the autocorrelation function K T ; b) the log-log 
plot of the normalized pdf of zero passage time of t f c 
supplemented by the local slope information. 



In this section we try to avoid the problem of finite- 
size-scaling related to true equilibrium critical temper- 
ature T c and critical exponents. The problem is solved 
by the multi-lattice self-organized simulations based on 
the dynamical rules treating the information from run- 
ning averages of magnetization. The considerations are 
addressed to models, where 2nd order phase transitions 
take place. The proposal is again applied to 2D Ising 
fcrromagnet on the square lattice. 

The quantity indicating deviations of magnetization 
order parameter tol = (l/£ 2 )X)(ij) s i lrom the gaus- 
sianity is the fourth-order cumulant 
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The standard way 

leading to true T c = limjv-»oo T C (N) is the construction 
of the temperature dependences Ul u t, Ul s ,t for two 
lattices L s ^ L\. Then T c follows from the condition of 
the scale-invariance 



Ul„t c = U Li 



(20) 



According it the self-organized multi-lattice MC sim- 
ulation method consists of the following three main 
points repeated in the canonical order for the counter 
*= 1,2,...: 



The performing of Lf spin flips on the lattice in- 
dexed by I and L 2 S spin flips on the second lat- 
tice. The flips are generated for fixed temperature 
T t . After it, the instant magnetizations (per site) 
m,L h T t and m,L 3 ,T t are calculated. 



2. The 



update 



of 



the cumulants (m 2 Li ),^ tiTt , (m 2 La } Vtt ,T t , (m 4 Ll }n,t.T t , 
{m\ ) v ,t.T t according to Eq.(ll), which yield to the 
modified of definition from Eq.(19) 
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3. The temperature shift 

T t+1 =T t +r t Asign(i$ a ) 
biased to eliminate difference 
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If Li > L s the ordering of cumulants in F^ ls is 
chosen subject to assumption 



U L „T t 
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T t >T c . 



(24) 



Any modification of the Eqs.(22) and (23) is possible 
when the preliminary recognition of the critical point 
neighborhood is performed. The parametric tuning 
recovers that stabilization of the noncanonical equi- 
librium via feedback F^ ls requires smaller -q and A 

than single-lattice simulations based on action of Ff . 
The Eq.(22) can be generalized for n r lattices, i.e. for 
n rp = ^(n r — 1) competing lattice pairs labeled by /, s: 

n A 
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n r X n r 
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(25) 



where l/n rp term rescales additive contributions. 

The presented method also offers the continuous 
checking of estimated critical exponents. It comes 
from the standard assumption that canonical equilib- 
rium magnetization exhibits critical scaling (\mL j \)t — 

LJ P/ "f m {Lj(T t - T c )), where / m (-) is the scaling func- 
tion and (3/v is the ratio of magnetization (/3) and the 
correlation length (v) critical exponents, respectively. 
If the temperature fluctuates nearly T c , the equilibrium 
finite-size-scaling relation changes to (\rriLj \) n ,t,T t =T c — 

Lj /m(0). For n T > 2 lattices and sufficiently small 
7], A, the following arithmetic average can be defined 

/ {\m Ll \)r,,t,T t \ 
V (|mL,|),,(,T, J 
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Similar to treatment of the steady temperature fluctu- 
ations, the quantity (3/v ~ N~^Y^t=i(P/ v )v-t can be 
defined. The simulations carried out for cases n r = 2, 3 
are compared in Fig. 4. In agreement with expectation, 
the localization of T c for n r = 3 with recursion taken 
from Eq.(25) is much subtle than for n r = 2. In ad- 
dition, the statistics of ((3/v) Vt t is weakly depending 
on rj and A, which seems to be logical due to univer- 
sality of exponents in the canonical equilibrium limit 
(T c = const). 

The n r = 2 simulations applied for L\ = 10, L 2 = 20, 
■q = 10~ 4 , A = 10~ 8 leads to the noncanonical equi- 
librium, for which the temperature average is associ- 
ated with estimate T c ~ 0.5667 of the exact value 
T c cx ~ [2 ln(l + V2]— 1 ~ 0.56729. The ratio (3/v ~ 0.122 
approximates the exact index (/3/^) cx = 0.125. Much 
slower walk for A = 10~ 9 provides only a insufficient 
improvement of the previous results. More appealing 
are estimates T c = 0.5673, (3/v ~ 0.123 obtained for 
n r = 3, Li = 10, L 2 = 20, L 3 = 30, rj = 10~ 4 , 
A = 10 -9 , A^av = 5 x 10 8 with balance of cumu- 
lants attained for ?7l=io,20,30,t c — 0.61. Note that 
[3/v does not change substantially [(3/v ~ 0.123(5)] 
if estimated from averages {\rriL 1 =w\)'q,t,T t = 0.37, 
(\m L2 =2o\) v ,t,T t = 0.34, < \m L3 =3o\)r,,t,T t = 0.33. 



5 The mean-field analysis of al- 
gorithm 

In this section we present calculations aimed to under- 
stand how the attractivity of critical point arises and 
how the noncanonical equilibrium is attained by means 




Figure 4: The stationary statistics of lattices coupled 
by the fourth-order magnetization cumulants. a) The 
pdf 's of temperature obtained for parameters rj = 10~ 3 , 
A = 10~ 4 . n r = 2 for L x = 10, L 2 = 20 [see (i)]. In 
that case additional low temperature bound T t > 0.1 is 
used to confine dynamics into the region T t > 0. For 
n T = 3, L\ = 10, L 2 = 16, L3 = 20 (ii) the stabilization 
bound is not necessary. For 77 = 10~ 4 and for A = 10~ 5 
the location of T c is much better [ see n r = 2 (hi) and 
n T = 3 (iv)]; b) The comparison of pdf's of the zero 
passage time t f u of two coupled lattice systems of sizes 
Li and L 2 - Calculated for 77 = 10~ 4 , A = 10~ 5 and 
T] = 10~ 4 , A = 10~ 8 . The log-log plot of pdf's results 
in the slope —0.671 in the region of vanishing t f u . For 
the middle Tpu region the slope is —0.59; c) pdf's of the 
effective critical index (/3/v) Vi t obtained for parametric 
choices from b). The arrow indicates the of exact T c . 



of feedback. Only a rough approximation of the com- 
plex simulation process is considered, where spin de- 
grees of freedom are replaced by the unique magnetiza- 
tion (per site) term m(t). Furthermore, it assumes that 
the selected central spin s(t) flips in a mean field created 
by its z neighbors. Let n(t) denotes the probability of 
the occurrence of s(t) = 1/2 state, then the probability 
of s(t) = —1/2 is 1 — n(t). The master equation for ir(t) 
can be written in the form 



J = (1 - n)W [ - +] - irW l+ - ] . 



(27) 



The Glauber's [24] heat bath dynamics with the tran- 
sition probabilities VF [_+1 and between states 
s(t) = ±1/2 is preferred in comparison to non- 
diffcrcntiable Metropolis form due to analyticity argu- 
ments relevant for formulation by means of differential 
equation. Within the mean-field approximation it can 
be assumed 



w [+ - ] (t) = 
w [ - +] {t) = 
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In the above expression Tf is the time associated with 
the spin flip process. The expression takes into account 
±zm(t) variations of energy belonging to flips from 
s = ±1/2 to s = ±1/2 given by the effective single-site 
Hamiltonian —zs(t) m(t). Assuming that it = m + 1/2 
and using Eqs.(27), (28) we obtain 
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Subsequently, the feedback differential equation of the 
temperature variable is suggested in the form 



dT 
dt 



a ( m 2 — m 2 ) = 



F m (t), 



(30) 



where m c is the "nucleation" parameter of the ferro- 
magnetic phase, and a > is the constant parameter. 
Unlike the works [18, 19], where feedback consisting of 
\m\ term is considered, the m 2 is absorbed to the feed- 
back F m {t) proposed here to ensure the analyticity. For 



to 2 > m 2 the temperature increases, whereas m 2 



< mt 



leads to the cooling. The stationary solution of Eqs.(29) 
and (30) is 



m c = ltanh(^) . 



(31) 



In the limit of vanishing m c , the solution of Eq.(31) can 
be written in terms of the inverse Taylor series in m c 
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where m c — corresponds to the known mean-field 



critical temperature T c 



MF 



z/4. Small negative shift 



of stationary T from Eq.(32) caused by m c ~ cor- 
responds to Fig. 5 including the numerical solution of 
Eqs.(29) and (30). 
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Figure 5: The numerical solution of differential equa- 
tions Eqs.(29) and (30) is presented. The transient dy- 
namics of the temperature obtained for m c — 0.005, 
a = 10~ 3 and initial conditions T| t=0 = 0.125, m\ t= o — 
0.05. 



6 The comparison 
SOC dynamics. 



of MC and 



In the section we discuss the universal aspects of the 
non-equilibrium self-organized MC dynamics applied to 
the canonical Ising model. As is already mentioned, 
the attributes of the SOC systems are avalanches re- 
flected by the power-law pdf distributions. We follow 
with the construction of certain temporal characteris- 
tics by supposing their uncertain links to avalanches. 
By using Eqs.(12)-(14) the evolution of any quantity 
can be mapped to the sequence of passage times. The 
example of this view represents pdf of Tpc depicted in 
Fig. 2. Because of the substantial difference in the expo- 
nents of pdf 's belonging to t f c and t f u , no universality 
attributes are indicated. More encouraging should be 
to find of pdf 's independent of the feedback type. The 
natural way toward this aim seems to be the investiga- 
tion of the passage time sequences linked to the order 
parameter of the canonical equilibrium of given system. 
In the case of the Ising model the ordering is described 
by the magnetization, or, eventually by isolated spin 
value. Therefore, it seems to be logical to define the pas- 
sage times T m , t s given by Eqs.(12)-(14) (corresponding 
to A = m,L,Si, where i is the arbitrary but fixed site 
position). The simulation results are depicted in Fig. 6. 
In their structure, the following attributes relevant for 
interpretation in terms of SOC can be identified: 

(I) the power-law behavior 

pdf(r s ) ~ , <f> ~ 1.3 



(33) 



with the unique exponent <f> pertaining to different 
feedbacks Ff , F^ , i.e. to single- lattice and two- 
lattice systems; 

(II) the interval of dependence from (I) broaden with 
the size of lattices 

(III) the exponent <j> (at the present level of accuracy) 
indistinguishable for pdf's taken for sequences 

{T^ k) } k =i and {rl fe) } fc= i. 
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Figure 6: The comparison of stationary pdf's of zero 
passage times t s and r m of different lattice sizes and 
different feedbacks. Simulated for r\ = 10 and A = 
10~ 5 . The figure shows rescaled pdf's of different sys- 
tems subject to settings ordered to three-component 
tuples: [feedback type, (quantity leading to {r^} se- 
quence), system size]: (i), Ff, (A — tul), L = 10: 
(ii), Ff , (A = Si ), L = 10; (iii), F t u , (A = m L ), L = 10; 
(iv), FY, (A = Si ), L = 10; (v), Ff, (A = Si ), L = 200; 
(vi), Ff, (A = St), L = 50; (vii), Ff , (A = 8i ), L = 4; 



The high-temperature limit of pdf of {ri fe ^}fc=i can be 
easily derived due to assumption about the absence of 
spin-spin correlations. Its form 



p- 



P=l/L 2 



(34) 



expresses invariance of s, during (r s — 1/L 2 ) spin flips, 
and its immediate change after the Tth spin flip occur- 
ring with probability p of the random picking of zth 
site. The simulations carried for paramagnet T 3> T c 
depicted in Fig. 7 agree with formula Eq.(34). 

Below T c pdf splits into separable contributions fit- 
ted here by the bimodal distribution P Ts = b V Ts (po) + 
b\P Ts (p\) with parameters b = 0.003, b\ = 0.004, 
Pa = 0.61, pi = 0.02. The supplementary analysis of 
statistics of the successive time differences of tol re- 
covers that biP Ts (pi) term originates from the mecha- 
nism of the long-time tunneling among nearly saturated 
states of the opposite polarity. From the figure it also 
follows that power-law short-time regime described by 
Eq. (33) is formed only if the feedback mechanism is ac- 
tivated. This conditional occurrence of universality can 
be considered as an additional (IV)th attribute relevant 
for identification of SOC. The noncanonical equilibrium 
attained by the self-organized MC dynamics for L = 10 
leads to P Ts dependence, which can be approximated 
by the fit 



5 2 r s v exp 



(35) 



with parameters b 2 = 0.00264, r 2 = 43.521. 

Let us to note that for the sand pile model [12] the 
spatial measure of avalanche is associated with the en- 
ergy integral taken during the stage following distur- 




Figure 7: The pdf of the spin passage times obtained 
for L = 10, f] = W- 4 , A = 10~ 6 . The MC simulations 
for fixed temperatures (i) T = 4T c (iV) (para-phase); 
black-dashed line corresponding to V Ts from Eq.(34); 
(iii) T = 0.8T c (iV) (ferro-phase) . The white-dashed 
line is the fit of the "bimodal" pdf V Te . Compared with 
self-organized MC simulations (ii), where feedback Ff 
yields to the power-law behavior with exponent from 
Eq.(33). 



bance. In the case of MC the analog of spatial measure 
can be the extremal magnetization 



'"max t = t (k) + 1 



max 



mL s ,T t 



(36) 



where t 



(k) 



The simulations show that pdf 

,(fc) 



corresponding to sequence {minLx}k=i can be approxi- 
mated by pdf(m max ) oc (u! max )^ m with (f> m = -2.0 ± 
0.1 obtained for the narrow span m max S< 0.2, 0.4 >. 
In agreement with SOC attribute labeled I, the inde- 
pendence of feedback type is identified. 

It should be noticed that there remains the principal 
problem of the link between SOC and presented self- 
organized MC method focused to the critical region. 
The problem is to identify the specific algorithmic seg- 
ment which should be interpreted as an analogue of 
a disturbance initializing the avalanche. Fortunately, 
the advanced MC approach exists through which a dis- 
turbance can be absorbed into MC algorithm with mi- 
nor violation of the original dynamics. In general, the 
approach of interest based on the coevolution of given 
system and its replica is known under the term damage 
spreading technique [25]. To apply it, let us consider 
the self-organized MC referential system labeled here 
as {l}t, which incorporates instant cumulants (single 
with Ff or multi-lattice based on Ff), spin configu- 
rations and instant temperature Consider also 
the replica counterpart {2} t of the system {1}*. As is 
known, the parallel simulation of {l} t and {2} t should 
be applied with the identical pseudorandom sequences. 
Canonically, the measure of damage effect is then de- 
fined through the single-time two-replica difference 



D 



{1,2} {1} 
L,,T, 



{1} 



l< 2 > 



{2} ) 



(37) 



where and m' 2 ' are magnetizations of two 

Li,T t {1} L u T t {2} 6 

lattices of the same size L\ belonging to systems {l} t 



and {2} t . Using Eqs.(12)-(14) the sequence of differ- 
t = 1,2,... is mapped onto the 



ences A* 



D 



{1,2} 



sequence of passage times Tp\ fc = 1,2, .... In the 



^2i=i T iP > tne re P nca {2jyo is rebuilded within two 
steps: 

1. all of the instant cumulants and spin configurations 
involved in {2} t are replaced by {l} t , i.e. D x ^ ) ' = 
0. 



l D 



2. the replica temperature is modified according 



T {2} _ T {i} 



(38) 



where constant €t causes the small disturbance 
of the coincidence of the contents of {l} + (fo and 

{2} t c*>. 

Evidently, the temperature disturbance plays role simi- 
lar to adding of the grain to the sand pile. The sand pile 
is stabilized if the rest state occurring for t — tp + 
is reached. The main idea of replica approach is that the 
relative motion of {2} t with respect to {l} t enhances 

-fl 2\ 

the nonlinearity responsible for a wide range of D] ' 

responses to ct- Thus, the only stochastic elements 

of replica simulation originate from the instants over 

which the content of {2} .Ah) is replaced by the content 
t—t D 

of {l} f _ f (fc). In analogy to Eq.(36), the complementary 
t—t D 

measure reflecting the spatial activity can be 



D (k) = 



D 



{1,2} 



(39) 



Using this, the simulated path is mapped onto the se- 
quence {-Dm2x}/c=i- Consequently, the pdf's can be ex- 
tracted which arc depicted in Fig. 8. They show that 
the effective exponents centred nearly —2.3 fit simulated 
pdf's fairly well. As in the case labeled I, pdf's related 
to rffl and Z?m2 x are weakly susceptible to the feedback 
choice. Since both temporal and spatial power-law at- 
tributes of universality are indicated, the standard SOC 
paradigm can be considered as a framework adaptable 
for the analysis of the suggested self-organized MC dy- 
namics. 



7 Conclusions 

Several versions of MC algorithm combining the self- 
organization principles with the MC simulations have 
been designed. The substantial feature of method is 
the establishing of running averages coinciding with 
gamma filtering of the noisy MC signal. The simu- 
lations are combined with the mean-field analysis de- 
scribing the motion of temperature near to magnetic 
transition point. The replica-based simulations indicate 
that pdf distributions of passage times in a noncanoni- 
cal equilibrium attain the interval of the power-law be- 
havior typical for the standard SOC pdf distributions. 
We hope that the present contribution will stimulate 




Figure 8: The stationary power-law pdf distributions of 
the passage times td calculated for different feedbacks 
(SOC criterion no. I) Ff, , self-organization param- 
eters r/ = 10~ 4 , A = 1CP 5 and disturbance parameter 
e T = 1CT 3 for sizes L l = 10, L 2 = 20. 



further self-organized studies of diverse lattice models, 
e.g. those related to percolation problem. 
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